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Abstract. Kinetic energies of a system of 4 He are investigated at zero temperature. 
The multi-weight extension to the diffusion Monte Carlo method is used to implement 
the Feynman-Hellmann theorem in an effective way. This method allows the quantities 
of interest to be computed with excellent accuracy. In order to study the importance 
of symmetry in the kinetic energy calculations, we have considered for the solid phase 
two guiding wave functions: the Nosanov-Jastrow without boson symmetry and the 
symmetric Nosanov-Jastrow with boson symmetry. In general very good agreement is 
found with the experimental data at both the liquid and solid phases. 
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In quantum systems single particle momentum distributions n(k) and their associated 
kinetic energies Ej~ are of great interest for both experimental and theoretical physicists. 
These quantities can characterize the extent in which these systems deviates from 
classical physics. The strong quantum effects present in the condensed phases of 4 He 
atoms offer an example where the n(k) can not be described by the Maxwell-Boltzmann 
distribution typical of classical systems. 

The most widely studied property of solid helium using neutrons is the kinetic 
energy [I]. To date, only experimental techniques associated with neutrons allow 
for direct measurements of the kinetic energy and the momentum distribution of this 
system. In the liquid phase these quantities have also been the aim of many theoretical 
studies. 

Estimations of the “exact” kinetic energy at zero temperature can be made 
through the diffusion Monte Carlo (DMC) method. However its computation is not 
straightforward as for the total energy. This is because the kinetic energy is a 
quantity that does not commute with the Hamiltonian of the system. A possible 
approach to estimate the kinetic energy is through the Hellmann-Feynman theorem 
using the multi-weight diffusion Monte Carlo method that allows the involved 

derivatives to be performed in a reliable way. In this extension of the standard DMC, 
the potential energy can be estimated without the usual difficulties associated with 
numerical differences of independent estimates and uncorrelated statistical uncertainties. 
In general, determinations of the kinetic energy might also be technically relevant in the 
context of projection Monte Carlo methods. This is because the total energy might have 
a bias proportional to the expectation value of the kinetic energy [4]. 

In this work we estimate the total and kinetic energies for the liquid and solid 
phase of a system formed by 4 He atoms. In the solid phase, we have explored in 
what degree two guiding functions with different overlaps with the true ground state 
would affect the results for the kinetic energy of the system. We have experienced with 
both a guiding function of the Nosanov-Jastrow (NJ) form and the so-called symmetric 
Nosanov-Jastrow (SNJ) function (q, h] ■ We want to verify how the introduction of a 
symmetric guiding function in the calculations would change the results obtained with 
a function that do not has this property. 

In the next section we describe the helium atoms system considered in this work. 
We also present the way in which the Feynman-Hcllmann theorem is implemented in 
the multi-weight diffusion Monte Carlo (DMC) method to compute kinetic energies. 
Results are given and some final comments are made in the last sections. 

2. Method 

We want to compute the kinetic energy of a system formed from 4 He atoms. With 
this aim we consider N bodies using periodic boundary conditions described by the 
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Hamiltonian 


H 
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N 
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where V (r^) is the interaction potential that depends on the inter-atomic distance 
r t j = |rj — Tj\. We have considered the hfd-b3-fcci1 potential proposed by Aziz and 
coworkers [7] 


V (x — r/r m ) 
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A e -( ax +Px 2 ) 
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3=0 
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where r m is the position of the potential minimum. The values of e, a, f3 and 6 * 21+6 can 
be found in [7j. The function B (x) is given by 


B(x) 


e if x < D, 

1 otherwise, 


(3) 


where D is another parameter of the inter-atomic potential. 

A standard method to determine the ground-state energy of this system is the 
DMC method M- This is accomplished by the simulation in imaginary time of a 
corresponding classical diffusion process with a source V(r). This process is guided by 
a trial wave function able to give a meaningful description of the ground state of the 
system. For efficiency, branching of configurations is introduced to avoid the burden of 
keeping those configurations that practically do not give any contribution to the results. 

In the solid phase we have experienced with two different guiding functions. One 
of them is of the Nosanov-Jastrow form 


^ NJ = ifj (R) <J> (R), (4) 

i’j (R) =II e 2 ^ rij) ’ ( 5 ) 

i<3 

N 

$(R) =Y[e-^ ri ~ ul2 , ( 6 ) 

2—1 

where R = {r l5 ■ • • ,rAr}, 1, are the lattice sites of the crystal, and finally b and c are 
variational parameters. A second guiding function was also considered in order to infer 
the importance of a symmetric guiding function in the calculations of the kinetic and 
potential energies. We choose to employ a symmetric Nosanov-Jastrow function given 
by 

t 57 VJ (r)=^(R) nE e " firi " ij|2 ’ 

3 i 

already used in DMC calculations 0 0 . It is important to note that this is an 
approximation to a true symmetric Nosanov-Jastrow wave function. In the liquid phase 
a Jastrow factor ifj( R) was used as a guiding function. 
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Since the kinetic energy operator does not commute with the Hamiltonian, its 
estimation is not straightforward in a DMC calculation. For this reason, we have 
implemented the Hellmann-Feynman theorem in the multi-weight extension of this 
method. The kinetic energy is then obtained from the estimated values of the potential 
and total energies. 

The potential energy E p is computed through the substitution V —> XV in the 
Hamiltonian and the total energy E x derivative with respect to the parameter A 


^ = s>< a »L- 

Finally the kinetic energy E k is computed through 


( 8 ) 


E k = E ( A=1 ) - E p . 


(9) 


In the implementation of the Hellman-Feynman theorem through the multi-weight 
extension to the DMC method m i, three different weights are associated to a given 
configuration or walker R, one for each value of A = {1 — 5, 1,1 + 5}. The weights 
are independently treated and the total energies of the systems are computed for each 
Hamiltonian H( A). Of course, a total energy corresponding to a particular value of A 
should agree within statistical uncertainties with the result obtained by the standard 
DMC simulation using the Hamiltonian H( A). 

In our calculations a single set of walkers is generated for all Hamiltonians H( A). 
The idea is to obtain the total energies E x with correlated statistical fluctuations. 
The values obtained in this way allow us to perform the numerical computation of 
the derivatives involved in the application of the Feynman-Hellmann theorem with the 
required accuracy. 

The drift of a given configuration R can be kept unique without any difficulty for 
all the values of A. A new configuration R' is sampled from R according the standard 
procedure through 


G d (R, R') 


( _I_^ 2 —^IR-R'-^AtDr ')! 2 

\4:7cDAt J 


( 10 ) 


where D = is the diffusion constant, At is the time step and F = 2Vln'F-r. 

In the calculation of the weights associated to a walker, required by the 
implementation of the Hellmann-Feynman theorem in the DMC method, a particular 
value of the parameter Ef is used for each value of A. Weights are updated by the factor 

Gl (R.R') = (11) 


when a configuration drifts from R to R'. In this expression Ef = H{X)^>t/^t are the 
local energies associated to the values of A we are considering. So, after the drift step, 
the weights are updated according to 


lu' x = oo x G x b (R, R'). 


( 12 ) 
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Table 1 . Liquid phase results for the total Ep, kinetic Ek and potential Ep energies 
in units of Kelvin at the given density. The guiding function was of the Jastrow form. 


p (^- 3 ) 

Et 

Ek 

E P 

0.3280 

-7.08±0.01 

11.96±0.03 

-19.04i0.03 

0.3423 

-7.20±0.01 

12.75i0.02 

-19.95±0.02 

0.3650 

-7.27i0.01 

14.44i0.04 

-21.71±0.04 

0.3874 

-7.21±0.01 

15.69±0.04 

-22.91±0.04 

0.4146 

-6.98±0.01 

17.75i0.04 

-24.74i0.04 

0.4380 

-6.57±0.02 

19.51±0.05 

-26.08±0.04 


The parameters Ef are periodically changed so that the sum of the configurations 
weights for that particular A is kept approximately constant. 

The energies at the £-th generation associated to each value of A are calculated in 
the usual way, as weighted averages over local energies of all the i walkers 


E} = 




E 


A 

i 


(13) 


Once the total energies are estimated, the kinetic and potential energies are computed 
through 

zpl+<5 rp 1 — 5 

V ‘= ‘ 25 * ’ K i = E t l -Vc (14) 

When the simulation reaches equilibrium these quantities are blocked and uncertainties 
estimated. In our calculations we have set S = 10“ 4 . 

So far all the additions we have performed to the DMC method are straightforward. 
Changes that require some elaboration concerns the combination of walkers, because 
those regarding either splitting or keeping their weights are very easy. We adopt the 
following rules, i) A walker is copied if all of its weights are greater than 2. Each copy 
bears half of the original weights. This step can be repeated if needed, ii) Walkers 
with at least one weight between a threshold w t hr and 2 are kept without changes, in) 
Walkers i and j with all weights smaller than Wthr are combined. For each weight A 
one of the usual rules of the literature is followed. With probability fff x the weight 

wf + Wj is attribute to the walker i and the value zero to the walker j. The sum of 
weights is attributed to the walker j with probability 1 — ^ A and the value zero to 

the walker i. For a particular A if the sum of the weights is equal to zero both walkers 
keep the value zero pTj. iv) Delete a walker if all its weights are equal to zero. 

The value of Wthr is chosen by analyzing a compromise between two conflicting 
requirements. Ideally Wthr needs to be small to avoid walkers with one or two of their 
weights equal to zero. These walkers are undesirable because they spoil the correlation 
we need to achieve for accurate numerical derivatives. On the other hand we do not want 
to carry walkers with too low weights along the simulation. These kind of walkers would 
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not contribute to the final result. In our simulations we have verified that Wthr = 0.3 
offers a good compromise between these two requirements. 

It is worthwhile mentioning the procedure adopted in the trial energies update. 
Since we have to deal with three values of Ef, it is useful to have an automatic scheme 
of updates. In this work the update was made every twenty generations using the 
following expression 


£t = Et + 


Co 


In 


T2 ~ Tl 


\ Wo )' 


(15) 


where the first term Ef is the total energy averaged over a block of the last twenty 
generations, Co is a parameter which smooths the fluctuations of the weights sum in 
the calculation of the trial energies, r 2 — r± is the elapsed “time” interval since the last 
update, W x is the sum of the cof weights of all walkers at “time” r 2 and W 0 is another 
constant. It is a target value for the sum of weights, the value we want to keep in the 
simulation at equilibrium. Its value is approximately equal to the number of walkers 
kept during the simulation. The update of the trial energies according to (15) was useful 
in the sense that we were able to perform runs where all the walkers did not have any of 
their weights equal to zero. The expression of (15) is a reminiscent from the estimation 
of the energy by the growth estimator, 


E X = E X + 


In 


(Ti) 


(16) 


T 2 -n \W x (r 2 ) / 

associated to the fluctuation of the weights sum in the elapsed “time” interval t 2 — T\. A 
“time” r has a simple relation with the performed number of generations M (attempts 
to move all particles) through r* = Mi At, where Ar is the time step used in (10). 


3. Results 

3.E Liquid phase 

Estimates of the total energy per 4 He atom are displayed in table [l] These results, as 
a function of the density, are shown in figure [T| as well. The curve in figure [l] is a third 
degree polynomial fit to the estimated energies in the variable (p — po)/po that depends 
on the fitted parameter p 0 , the equilibrium density. Our fit gives p 0 = 0.366 ± 0.002, 
in excellent agreement with experiment. In general the total energies are also in good 
agreement with experimental data obtained at finite temperatures. A possible bias in the 
results (that produces an underestimation of the experimental data) can be attributed 
to the lack of three-body interactions in the inter-atomic potential we have considered 
in our calculations U21 US]- In figure [2] we show a comparison between the inter-atomic 
potential used in this work and the one dubbed hfdhe2 ra widely used in the literature. 
The former inter-atomic potential is more attractive than the last one and certainly a 
difference in energy can be attributed to the characteristics of these interactions. 

As we have discussed the total, potential and kinetic energies at a given density were 
estimated using correlated estimates of the desired quantities. In figure [3] we present 
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Figure 1. Equation of state as a function of the density in the liquid phase. 
Theoretical results are displayed by • . The curve is a fit to the estimated values 
of the energies (see text). The errors are smaller than the symbol size. Experimental 
data for the liquid phase [13 US] are shown by □ and A symbols. 



r (nm) 


Figure 2. Comparison of the inter-atomic potential hfd-b3-fci1 [71 employed in 
this work (dubbed Aziz 95 at the Fig. legend) with the widely used HFDHE2 [0] (Aziz 
79 at the Fig. key). 


the kinetic energy results of table [TJ The estimated values were fitted to a parabola. 
We considered a quadratic fit because this behavior was found in experimental data 
obtained at finite temperature czt Even if we have in mind that this behavior was 
determined above the A-transition. The experimental values displayed in the same 
figure were determined at T=0.05 K [18] . 

3.2. Solid phase 

A defect free hep crystalline structure with N = 288 4 He atoms simulation cell was used 
to estimate the total energy per atom in the solid phase. In figure [4] the results (as a 
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Density (o' 3 ) 

Figure 3. Kinetic energy per atom as a function of the density in the liquid phase. 
The curve stands for a fit to the estimated values displayed as •. Errors are smaller 
than the symbols size. The symbols □ represent experimental results [15] . 


function of the density) are shown by curves fitted to the values of the total energies 
displayed in table[2j In this phase, at the level of precision we have considered, the results 
obtained with a guiding function of the Nosanov-Jastrow form are indistinguishable 
within statistical uncertainty or are marginally lower than those determined with the 
symmetric Nosanov-Jastrow guiding function. The only exception is at p = 0.5277cr~ 3 
where the SNJ guiding function gives a marginally lower energy. The overall situation 
might be attributed to the different degrees of superposition that the guiding functions 
have with respect to the true ground state of the system. The equation of state is also 
displayed by curves that are third degree polynomial fits to the estimated energies in 
the variable (p — po)/po that depends on the fitted parameter p 0 . The parameter p 0 
does not have any particular meaning in this phase. Once again the total energies are in 
general in good agreement with the experimental data obtained at finite temperatures. 
The same bias observed in the liquid phase can be seen in the solid phase, and here 
also it is attributed to the lack of three-body interactions in the inter-atomic potential 
considered H21H3I- 

Correlated estimates of the quantities needed to estimate the kinetic energy were 
used and the values obtained are presented in figure [5j In this phase a straight line 
was used to fit the kinetic energy results of table [2] We considered a linear fit since 
this behavior was found in experimental data obtained at temperatures above the A- 
transition S3- The experimental values displayed in the figure are from mmm- 
Although in the solid phase the kinetic energies obtained with a Nosanov-Jastrow 
guiding function could be indistinguishable (within statistical uncertainty) from those 
of a symmetric Nosanov-Jastrow at a few densities, always the results from the first 
guiding function were lower than those obtained with the second one. This trend is not 
obvious, since atoms are in principle much less localized when the simulation is guided 
by the symmetric Nosanov-Jastrow guiding function. 
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Table 2. Total Ep, kinetic Ek and potential Ep energies in units of Kelvin for the 
solid phase at the given densities. The results were obtained with the Nosanov-Jastrow 
(NJ) and the symmetric Nosanov-Jastrow (SNJ) guiding functions. 


p (c Ep Ek Ep 



NJ 

SNJ 

NJ 

SNJ 

NJ 

SNJ 

0.5026 

-5.74±0.01 

-5.70±0.01 

26.37±0.04 

26.61±0.05 

-32.11d=0.04 

-32.31±0.05 

0.5126 

-5.54±0.01 

-5.53±0.02 

27.28±0.04 

27.40±0.10 

-32.82±0.03 

-33.00±0.10 

0.5277 

-5.16±0.01 

-5.20±0.01 

28.66±0.06 

28.73±0.05 

-33.83±0.06 

-33.93±0.04 

0.5344 

-4.98±0.01 

-4.93±0.01 

29.37±0.03 

29.85±0.05 

-34.36±0.03 

-34.78zt0.05 

0.5494 

-4.51±0.01 

-4.50±0.02 

30.76±0.07 

30.99±0.07 

-35.28±0.08 

-35.49±0.06 

0.5694 

-3.78±0.03 

-3.70±0.01 

32.80±0.10 

33.15±0.05 

-36.60±0.10 

-36.86±0.05 

0.5895 

-2.85±0.01 

-2.80±0.01 

34.81±0.02 

35.10±0.10 

-37.67±0.02 

-37.90±0.10 



Density (o' 3 ) 

Figure 4. Total energy per atom as a function of the density in the solid phase. 
Theoretical results are displayed by curves fitted to the estimated values of the 
energies (see text). Results obtained with the NJ (blue line) and the SNJ (dashed 
and thick green line) guiding functions are hard to distinguish at the figure resolution. 
Experimental data mmm are displayed by O, their errors have been assumed to 
be half of the penultimate less significative figure. 


4. Final Comments 

The multi-weight DMC method has shown to be reliable to compute the kinetic energies 
of systems formed by helium atoms. Furthermore it is not hard to realize that there 
are other applications of the multi-weight DMC method. One of them in qhich we are 
working is its application in the refinement of guiding functions of systems that obey 
Fermi statistics without the need of variational approaches. This would be possible if 
small changes in a guiding function could be mapped into changes of the inter-atomic 
potential. We already know that the multi-weight method can be used to analyze the 
consequences of small changes in the interacting potential mm- 

In the liquid phase as expected the agreement with experiment is very good. The 
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Density (o' 3 ) 

Figure 5. Kinetic energy per atom as a function of the density in the solid 
phase. The curves stand for fits to the estimated values (see text). Experimental 
data nu unusii na are displayed by O. Fits to the estimated values obtained with the 
NJ (blue line) and the SNJ (dashed and thicker green line) guiding functions almost 
coincide. 

kinetic energy in this phase at all densities we have considered are also in good agreement 
with experiment and other theoretical estimates of this quantity. 

Our results in the solid phase show that both guiding functions, the Nosanov- 
Jastrow and the symmetric version [6] we have considered, can give energies in good 
agreement with experiment. However, the trend we have observed for the estimated 
values of the energies suggests that the Nosanov-Jastrow guide function can be the 
preferred form when a simple DMC calculation needs to be performed. Although the 
differences in the results obtained with these guiding functions might be attributed to 
the particular and well know fact that exchange is not so important in solid helium, 
the result could also be a consequence of the approximation used in the symmetric 
Nosanov-Jastrow function. It does not have all the terms that a true symmetrization 
would give. 

In conclusion, the extension of the DMC method we have applied in this work 
allows an accurate determination of quantities that depends on differences of estimates 
even when these differences are very small. The method avoids the use of quantities 
with uncorrelated statistical uncertainties to compute their differences. This is a well 
known difficulty of Monte Carlo methods that can lead to a complete loss of accuracy of 
the quantities of interest. Moreover and more important, the application of this method 
in slightly different contexts may bring an enhancement of our understanding of the 
quantum many-body systems. 
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